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SMALL SCALE CLUSTERING IN THE ISOTROPIC ARRIVAL DISTRIBUTION OF ULTRA-HIGH 
ENERGY COSMIC RAYS AND IMPLICATIONS FOR THEIR SOURCE CANDIDATES 
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ABSTRACT 

We present numerical simulations on the propagation of UHE protons with energies of (10 9 - 



10 22 ) eV in 



extragalactic magnetic fields over 1 Gpc. We use the ORS galaxy sample, which allow us to accurately quantify 
the contribution of nearby sources to the energy spectrum and the arrival distribution, as a source model. The 
sample is corrected taking the selection effect and absence of galaxies in the zone of avoidance (\b\ < 20°) into 
account. We calculate three observable quantities, cosmic ray spectrum, harmonic amplitude, and two point 
correlation function from our data of numerical simulations. With these quantities, we compare the results of our 
numerical calculations with the observation. We find that the arrival distribution of UHECRs become to be most 
isotropic as restricting sources to luminous galaxies (Mn m = -20.5). However, it is not isotropic enough to be 
consistent with the AGASA observation, even for Mi; m = -20.5. In order to obtain sufficiently isotropic arrival 
distribution, we randomly select sources, which contribute to the observed cosmic ray flux, from the ORS sample 
more luminous than -20.5 mag, and investigate dependence of the results on their number. We show that the 
three observable quantities including the GZK cutoff of the energy spectrum can be reproduced in the case that 
the number fraction ~ 10 _1 7 of the ORS galaxies more luminous than -20.5 mag is selected as UHECR sources. 
In terms of the source number density, this constraint corresponds to ~ 10~ 6 Mpc" 3 . However, since mean number 
of sources within the GZK sphere is only ~ 0.5 in this case, the AGASA 8 events above 10 20 eV, which do not 
constitute the clustered events with each other, can not be reproduced. On the other hand, if the cosmic ray flux 
measured by the HiRes, which is consistent with the GZK cutoff, is correct and observational features about the 
arrival distribution of UHECRs are same as the AGASA, our source model can explain both the arrival distribution 
and the flux at the same time. Thus, we conclude that large fraction of the AGASA 8 events above 10 20 eV might 
originate in the topdown scenarios, or that the cosmic ray flux measured by the HiRes experiment might be better. 
We also discuss the origin of UHECRs below 10 20() eV through comparisons between the number density of 
astrophysical source candidates and our result (~ 10" 6 Mpc" 3 ). 



Subject headings: cosmic rays — methods: numerical - 
large-scale structure of universe 

1. INTRODUCTION 

In spite of a increasing amount of data, the origin of cos- 
mic rays in particular at the highest energy (~ 10 20 eV) is still 
unknown. It is very difficult to produce such ultra-high en- 
ergy cosmic rays (hereafter UHECRs) in astrophysical shocks, 
which are thought to be responsible for the galactic cosmic rays. 
The gyroradii of UHECRs with energies of 10 20 eV are larger 
than our Galaxy, and this suggests an extragalactic origin for 
any astrophysical scenario. Using a Hillas-plot (Hillas 1984; 
Selvon 2000), active galactic nuclei (AGN), gamma-ray bursts 
(GRBs), dead quasars and colliding galaxies are considered as 
probable candidates. These conventional acceleration models 
are called bottom-up scenarios. On the other hand, this situa- 
tion has in recent years triggered a number of exotic production 
models based on new physics beyond the Standard model of 
the particle physics (see Bhattacharjee and Sigl 2000, and ref- 
erences therein), which are called top-down scenarios. In such 
scenarios, UHECRs owe their origin to decay of some super- 
massive particles, which could be produced from Topological 
Defects (TDs), or be certain metastable supermassive relic par- 
ticles (MSRPs) produced in the early universe. 

In bottom-up scenarios in which protons are accelerated in 
the potential UHECRs sources, like GRBs and AGNs, the 
large distances to the earth lead to another problem, called 
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the GZK effect (Greisen 1966; Zatsepin and Kuz'min 1966). 
UHE protons with energy above 8 x 10 19 eV interact with pho- 
tons of Cosmic Microwave Background (CMB) and lose their 
energy rapidly due to the production of pions over distances 
of tens of Mpc (Stanev et al. 2000). This should results in a 
cut-off of the cosmic ray spectrum at energy around 8xl0 19 
eV. The situation for nuclei is considered to be worse due to 
the photo-disintegration mechanism (Stanev et al. 2000). Ob- 
servationally, there seems to be a disagreement between the 
AGASA which revealed the extension of the cosmic ray en- 
ergy spectrum above 10 20 eV (Takeda et al. 1998) and the High 
Resolution Fly's Eye (HiRes; Wilkinson et al. 1999) which 
recently reported the cosmic ray flux with the GZK cut-off 
(Abu-Zayyad et al. 2002). 

Since UHECRs above 10 20 eV observed by the AGASA 
must originate within the GZK sphere, whose radius is typi- 
cally about 50 Mpc, the arrival directions are expected to point 
toward their sources within a few degrees, which is typical de- 
flection angle of UHECRs (e.g., Blasi and Olinto 1998). How- 
ever, no plausible source counterparts have been found within 
the GZK sphere. There is also a problem concerning the ar- 
rival directions of UHECRs (e.g., Blasi and Olinto 1998). Their 
distribution seems to be isotropic on a large scale with a sta- 
tistically significant small scale clustering (Takeda et al. 1999). 
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The current AGASA data set of 57 events above 4 x 10 eV 
contains four doublets and one triplet within a separation angle 
of 2.5°. Chance probability to observe such clusters under an 
isotropic distribution is about 1 % (Hayashida et al. 2000). Po- 
tential models of UHECR origin are constrained by their ability 
to reproduce the observed energy spectrum and the isotropic ar- 
rival distribution of UHECRs with the small scale clustering. 

A number of studies has been devoted to the problem of the 
energy spectrum. Yoshida and Teshima (1993) show the GZK 
cutoff in consideration of evolution of the universe. Blanton 
et al.(2001) show that the observed number of events above 
10 20 eV is 2 a higher than the expected one, assuming a hard 
(E~ 2 ) injection spectrum and small local overdensity (a factor 
of two). Berezinsky et al.(2002) analyze three source models, 
GRB, AGN and local enhancement of UHECR sources. They 
conclude the AGASA events above 10 20 eV, if confirmed, must 
be of different origin. Although these authors show that it is dif- 
ficult to reproduce the observed extension of the energy spec- 
trum by the bottom-up scenarios, they do not take into account 
of the effects of the extragalactic magnetic field (EGMF) on 
the energy spectrum. Furthermore, since they do not consider 
angular source distribution as well as the EGMF, the problems 
concerning the arrival distribution of UHECRs can not be dis- 
cussed. 

Although arrival distributions from a single source are cal- 
culated in the presence of magnetic field in some papers 
(Lemoine et al. 1997; Stanev et al. 2000), it appears to be dif- 
ficult to explain the observational isotropy by a single source. 
With consideration of both the magnetic field and the source 
distribution, the energy spectrum and the distribution of the ar- 
rival directions are calculated in Sigl et al.(1999) and Lemoine 
et al.(1999). However the source distribution seems to be too 
simple, a two dimensional sheet (Sigl et al. 1999) or a pancake- 
like profile (Lemoine et al. 1999), which represent the Local 
Super Cluster (LSC). 

Smialkowski et al.(2002) calculate arrival distribution of 
UHECRs using the IRAS PSCz catalog of IR galaxies as a 
realistic source model. However, they restrict their attention 
to galaxies with far infrared luminosity L F i R > 10 n L Q trig- 
gered by collision and merging processes, which are possible 
sites of UHECR acceleration. It is also suspected that this 
galaxy sample do not contain some nearby galaxies in the LSC. 
Since nearby sources within the GZK sphere are responsible for 
UHECRs above 10 20 eV, galaxy sample with better complete- 
ness on nearby galaxies should be used. 

Given these situations, we examined in our previous paper 
(Ide et al. 2001) whether the current observation can be ex- 
plained by a bottom-up scenario in which source distribution 
of UHECRs is proportional to that of galaxies. We used the 
Optical Redshift Survey (ORS; Santiago et al. 1995) as a re- 
alistic source model. Completeness on nearby galaxies of this 
sample would be better than that of the IRAS PSCz. By using 
the ORS sample, we can accurately quantify the contribution 
of nearby sources to the energy spectrum above 10 20 eV. We 
also took effects of the EGMF on the energy spectrum and the 
arrival distribution of UHECRs into account. 

However, there are problems in our previous study. First, we 
did not correct absence of galaxies in the zone of avoidance 
(\b\ < 20°) when using the ORS sample, so that we could not 
present accurate discussion of statistics on the arrival directions 
of UHECRs. Second, we considered propagation of UHECRs 
only over 40 Mpc, focusing on UHECRs with energies above 



8 x 10 eV. In this case, contribution of distant sources (out- 
side 40 Mpc) to the cosmic ray flux below 8 x 10 19 eV can not 
be taken into account. Finally, we injected UHECRs from all 
the galaxies including dwarf ones. Luminous galaxies are gen- 
erally larger than dwarf ones, and thus they are expected to be 
able to accelerate particles to higher energy. Thus, we should 
investigate dependence of the results on the limiting magnitude 
of galaxies. 

This paper is an extended study of our earlier work. We cal- 
culate the energy spectrum and the arrival distribution of UHE- 
CRs using the ORS galaxy sample. We correct the sample tak- 
ing the selection effect and absence of galaxies in the zone of 
avoidance (\b\ < 20°) into account. This allow us to accurately 
discuss statistics on the arrival directions of UHECRs. We also 
extend the propagation distance of UHECRs, from 40 Mpc to 
1 Gpc. As s result, we can calculate the cosmic ray flux be- 
low 8 x 10 19 eV, properly including contribution from distant 
sources. Furthermore, we investigate dependence of the results 
on the limiting magnitude of galaxies. 

We compare the results of the energy spectra with the cur- 
rent observations, AGASA and HiRes. In order to explore the 
possibility of reproducing the extension of energy spectrum be- 
yond 10 20 eV, we adopt several strength of the EGMF (1,10 
and 100 nG). We also investigate dependence of the arrival dis- 
tribution of UHECRs on the limiting magnitude of galaxies and 
the source number, and estimate the most favorable ones from 
comparison with the AGASA data using the harmonic analy- 
sis and the two point correlation function as statistical meth- 
ods. We show that the three observable quantities including 
the GZK cutoff of the energy spectrum can be reproduced in 
the case that the number fraction ~ 10 _1 7 of the ORS galaxies 
more luminous than -20.5 mag is selected as UHECR sources. 
Implications of the results for the UHECR source candidate are 
discussed in detail. 

In section 2, we show our method of calculation. Results 
are shown in section 3 and their implications are discussed in 
section 4. We conclude in section 5. 

2. METHOD OF CALCULATION 
2. 1 . Method of Calculation for Propagation of UHECRs 

This subsection provides the method of Monte Carlo simula- 
tions for propagating protons in intergalactic space. At first, we 
assume that the composition of UHECRs is proton, and inject 
an E~ 2 spectrum within the range of (10 19 5 - 10 22 )eV. 10000 
protons are injected in each of 26 energy bins, that is, 10 bins 
per decade of energy for the cases of l c = 40Mpc and (B, l c ) = 
(l,10nG, lOMpc), where B and l c are strength and correlation 
length of the extra galactic magnetic field (explained below). 
For another cases, we inject 5000 protons in each of the bins. 

Particles below ~8x 10 19 eV lose their energies mainly 
by pair creations and above it by photopion production 
(Yoshida and Teshima 1993) in collisions with photons of the 
CMB. The pair production can be treated as a continuous loss 
process considering its small inelasticity (~ 10~ 3 ). We adopt 
the analytical fit functions given by Chododowski et al. (1992) 
to calculate the energy loss rate for the pair production on 
isotropic photons. According to them, the energy loss rate of 
a relativistic nucleus for the pair production on isotropic pho- 
tons is given by 

--77 = ar o cZ — / dnn[ — \ — r , (1) 
dt m A J 2 \2~fJ k 2 
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where 7 is the Lorenz factor of the particle, k = 2kj ( k is the 
momentum of the particle in units of m e c), n(n), a, rrj, Z, and 
rtiA are the photon distribution in the momentum space, the fine- 
structure constant, the classical electron radius, the charge of 
the particle, and the rest mass of the particle, respectively. For 
the energy range E > 10 19 eV, <p can be well represented as 

3 

ip — > K^ffjln'K, (2) 
1=0 

d Q ~ -86.07, di ~ 50.96, d 2 ~ -14.45, d 3 = 8/3. (3) 

Since protons lose a large fraction of their energy in the pho- 
topion production, its treatment is very important. We use the 
interaction length and the energy distribution of final protons 
as a function of initial proton energy which is calculated by 
simulating the photopion production with the event generator 
SOPHIA (Mucke et al. 2000). We checked that arrival spec- 
trum at Earth for mono-energetic injection of energy E = 10 21 5 
eV is consistent with that of Stanev et al.(2000) (Fig. 2), which 
is also calculated with the event generator SOPHIA, for the case 
of (B, l c ) = (InG, IMpc) which is the closest parameter set to 
the one of Stanev et al.(2000). 

The EGMF are little known theoretically and observation- 
ally. The upper limit for its strength and correlation length 

fiiGM^y 2 < lnG(lMpc) 1 / 2 , as measured by Faraday rotation of 
radio signals from distant quasars (Kronberg 1994), is often 
used (Stanev et al. 2000; Smialkowski et al. 2002). We adopt 
in this study not only this widely accepted value, InG, but also 
relatively strong magnetic field (lOnG, lOOnG), which is bet- 
ter for explaining the extension of the energy spectrum and the 
observed isotropic arrival distribution. In this case, deflection 
angles of UHECRs become so large that the no counterparts 
problem can be simply solved. 

We assume that the magnetic field is represented as the Gaus- 
sian random field with zero mean and a power-law spectrum. 
Thus, (B 2 {k)) can be written as 

(B 2 (k)) cck" H for 2tt// c < A; < 27r// cut (4) 
= otherwise, (5) 

where l c is the correlation length of the EGMF and / cut charac- 
terizes the numerical cut-off scale. We use n H = -1 1 /3 corre- 
sponding to the Kolmogorov spectrum. According to our pre- 
vious study (Ide et al. 2001), the correlation length is chosen to 
be 1 Mpc, 10 Mpc, and 40 Mpc. 1 Mpc is roughly equal to mean 
separation of galaxies and widely used as the typical value of 
the correlation length. As we show later, the small scale clus- 
tering observed by the AGASA is difficult to be reproduced for 
l c > 10 Mpc (See Figure 13 and 17). We do not have to con- 
sider the correlation length larger than 40 Mpc. Physically one 
expect l cul <C l c , but we set / cut = 1/8 x / c in order to save the 
CPU time. 

The universe is covered with cubes of side l c . For each of 
the cubes, Fourier components of the EGMF are dialed on a 
cubic cell in wave number space, whose side is 2w/l c , with ran- 
dom phases according to the Kolmogorov spectrum, and then 
Fourier transformed onto the corresponding cubic cell in real 
space. We create the EGMF of 20 x 20 x 20 cubes of side l c , 
and outside it, adopt the periodic boundary condition in order 
to reduce storage data for magnetic field components. Simi- 
lar methods for the turbulent magnetic fields have been adopted 
(Sigl et al. 1999; Lemoine et al. 1999; Isola et al. 2002). 



In this study, we assume that the source distribution of UHE- 
CRs is proportional to that of the galaxies. We use the re- 
alistic data from the ORS (Santiago et al. 1995) galaxy cata- 
log, which is nearly full sky survey, and contains two subcat- 
alogs, one complete to a B magnitude of 14.5 (ORS-m) and 
the other complete to a B major axis diameter of 1/9 (ORS- 
d). Another nearly full sky catalog which is often used in this 
kind of study is the IRAS PSCz Survey (Blantonet al. 2001; 
Smialkowski et al. 2002). However, it is suspected that some 
galaxies in the LSC are excluded. Since local overdensity of 
UHECR sources makes the GZK cutoff less sharp or eliminate 
it, we use the ORS as a source model rather than the IRAS. 

In any statistical acceleration mechanism, particles must be 
kept confined within the acceleration site. Luminous galaxies 
are generally larger than dwarf ones, and thus they are expected 
to be able to accelerate particles to higher energy. Furthermore, 
it is less well known that luminous galaxies in the LSC dis- 
tribute outward than faint galaxies, contrary to general clus- 
ters of galaxies (Yoshiguchi et al. 2002). Accordingly, we in- 
ject UHECRs from galaxies with different limiting magnitudes, 
and explore the source model which reproduces the current ob- 
servation. As an example, this selection of luminous galax- 
ies corresponds to restricting the sources of UHECRs to host 
galaxies of AGNs. It is still unknown how much an ultimate 
UHECR source contribute to the observed cosmic ray flux. We 
thus consider the two cases in which all galaxies are the same, 
and they inject cosmic rays proportional to their absolute lumi- 
nosity. This allow us to investigate effects of the uncertainty of 
the cosmic ray luminosity on the results. 
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FIG. 2. — Absolute magnitude plotted as a function of the recession velocity 
for galaxies in our ORS sample. 

Since we consider the galaxies with different absolute mag- 
nitudes, we use the magnitude-limited subcatalog, i.e., ORS- 
m, which is also composed of two subsamples, ESO-m (2437 
galaxies) and UGC-m (3279 galaxies). We further restrict our- 
selves to galaxies within 8000 km s" 1 and with Mb -51ogio/? < 
-0.5-51ogiov, which corresponds to the apparent magnitude 
limit of m B = 14.5, where h is the dimensionless Hubble con- 
stant h = H()/( 100 km s" 1 Mpc" 1 ) and vis the recession velocity. 
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FIG. 1 .— Distribution of galaxies in the ESO-m and the UGC-m samples before corrected (v < 8000 km s" 1 , M B < -14.0, m B < 14.5, \b\ > 20°). 



In the following, we abbreviate M B - 51og ioh to M B . We exclude 
very faint nearby galaxies with M B > -14.0. Galaxies in the 
zone of avoidance {\b\ < 20°) are not observed. Our final ORS 
sample contains 5178 galaxies with 2346 galaxies in the ESO- 
m and 2832 galaxies in the UGC-m. Figure 1 shows the sky dis- 
tribution of our ORS galaxies in the equatorial coordinate. The 
strip between the ESO and UGC region (-17.5° < 5 < -2.5°) 
is covered by the ESGC, which is a part of the diameter-limited 
subcatalog, ORS-d. The absolute magnitudes are shown in Fig- 
ure 2 as a function of the recession velocity. Throughout the 
paper, we assume that there are no departures from uniform 
Hubble expansion for simplicity, and use h = 0.75 in calculat- 
ing distance between galaxies in the sample and our Galaxy 
from recession velocity. 

In order to calculate the energy spectrum and the distribu- 
tion of arrival directions of UHECRs realistically, there are 
two key elements of the galaxy sample to be corrected. First, 
galaxies in a given magnitude-limited sample are biased trac- 
ers of matter distribution because of the flux limit. For each 
galaxy in the sample, we compensate galaxies which is not 
included in the sample using the selection function given in 
Santiago et al.(1996). The positions of compensated galax- 
ies are determined according to Gaussian distribution whose 
mean is the position of the original galaxy and rms is such 
that a = 1/3 x (n^(r)) -1 / 3 , where r is the distance from our 
Galaxy and n, <p(r) are number density of galaxies which is also 
calculated using the selection function, value of the selection 
function at distance r from us respectively. For the reason of 
the selection effect, the ORS is only sampling the universe out 
to 8000 km s" 1 (see Figure 2). We assume that contributions 
from sources outside it is completely isotropic, and calculate 
its amount from the number of galaxies inside it. Second, our 
ORS sample does not include galaxies in the zone of avoidance 
(\b\ < 20°) and the ESGC region (-17.5° < 6 < -2.5°). In the 
same way, we assume that the source distribution in this re- 
gion is homogeneous, and calculate its number density from the 
number of galaxies in the observed region, which is corrected 
for the selection effect. In Figure 3, we show the sky distri- 



bution of the ORS galaxies corrected in the manner explained 
above. 




FIG. 3. — Distribution of galaxies in our ORS sample within 8000 km s , 
which are corrected for the selection effect in the manner explained in the text. 

Finally, we now explain how the angular probability distri- 
bution of UHECRs are calculated. At first, protons with a flat 
energy spectrum are injected isotropically at a given point in 
each of the 26 energy bins, and then propagated in the EGMF 
over 1 Gpc for 15 Gyr. Weighted with a factor corresponding 
to a E~ 2 power law spectrum, this provides distribution of en- 
ergy, deflection angle, and time delay of UHECRs as a function 
of the distance from the initial point. With this distribution, 
we can calculate the flux and the arrival directions of UHECRs 
injected at a single UHECR source. Then, summing contribu- 
tions from all the sources, we obtain the angular probability 
distributions of UHECRs. In this paper, we use the distribu- 
tion of energies and deflection angles integrated over the time 
delay, assuming that the cosmic ray flux at the earth is station- 
ary. In Figure 4, we show number of protons with energies of 
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(10 196 - 10 20 ' 3 ) eV integrated over deflection angle, time delay 
as a function of propagation distance. It is noted that this num- 
ber is proportional to the cosmic ray flux in this energy range 
from sources in concentric shells with an infinitely small width, 
assuming that a source distribution is homogeneous. 
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FIG. 4. — Number of protons with arrival energies of (10"- s - 10 203 ) eV in- 
tegrated over deflection angle and time delay after propagation from an initial 
point as a function of the propagation distance. The number is normalized to 
be 1 at the injection. It is noted that this number is proportional to the cosmic 
ray flux in this energy range from sources in concentric shells with an infinitely 
small width, assuming that a source distribution is homogeneous. 



2.2. Statistical Methods 

In this subsection, we explain the two statistical quan- 
tities, the harmonics analysis for large scale anisotropy 
(Hayashida et al. 1999) and the two point correlation function 
for small scale anisotropy. The harmonic analysis to the right 
ascension distribution of events is the conventional method to 
search for global anisotropy of cosmic ray arrival distribution. 
For a ground-based detector like the AGASA, the almost uni- 
form observation in right ascension is expected. The m-th har- 
monic amplitude r is determined by fitting the distribution to a 
sine wave with period 2ir/m. For a sample of n measurements 
of phase, <j>\, (f>2, ■ ■ ■, (f>„ (0 < 0; < 2ir), it is expressed as 

r = (a 2 + b 2 ) 1/2 (6) 
where, a = . cosmdi, b = , sinm</> ; . We calculate the 
first and second amplitude from a set of events generated ac- 
cording to predicted probability density distribution of arrival 
directions of UHECRs. 

If events with total number n are uniformly distributed in 
right ascension, the chance probability of observing the ampli- 
tude > r is given by, 

P = exp(-k), (7) 

where 

k = nr 2 /4. (8) 



The AGASA 57 events is consistent with isotropic source dis- 
tribution within 90 % confidence level as shown later. Accord- 
ingly, we adopt the amplitude such that P = 0. 1 as observational 
constraint, rather than the most probable amplitude calculated 
from the observed 57 events. 



The two point correlation function N(9) contains information 
on the small scale anisotropy. We start from a set of generated 
events or the actual AGASA data. For each event, we divide 
the sphere into concentric bins of angular size A9, and count 
the number of events falling into each bin. We then divide it by 
the solid angle of the corresponding bin, that is, 



N(9) = 



1 



2tt|cosi9-cos((9+A(9)| 



]T 1 [sr" 1 ], (9) 



9<4><0+A8 



where tfi denotes the separation angle of the two events. A9 is 
taken to be 1° in this analysis. The AGASA data shows strong 
correlation at small angle (<~ 2°) with 5 a significance of devi- 
ation from an isotropic distribution as shown later. 

In order to quantify the statistical significance of deviations 
between models and data we introduce x flli , which is defined 

by, 



l 



X^max. 



A 



ft=0 



(N(9 1 )-N ohs (9 i )) z 



(10) 



where -/V bs(#i) is the value of the two point correlation func- 
tion obtained by the AGASA observation at angle 9[, and a x is 
statistical error of the N(9\) due to the finite number of simu- 
lated events, which is set to be the observed one. In this study, 
we take 9 m . dx to be 10 deg for properly quantifying deviations 
of models from the sudden increase of correlation at the small 
angle scale observed by the AGASA. 

We compare the model predictions with the existing AGASA 
data throughout the paper. Therefore, we set the number 
of simulated events to be that of the AGASA data above 
10 19 - 6 eV (57), and restrict the arrival directions of UHE- 
CRs in the range of -10° < S < 80° when calculating these 
statistical quantities. There are also the present working 
or development of large-aperture new detectors, such as the 
HiRes (Wilkinson et al. 1999) and South and North Auger 
(Capelle et al. 1998). If they report the arrival directions of ob- 
served events, we can make comparisons between the model 
predictions and their data using the same statistical methods. 

3. RESULTS 
3.1. Angular Images of UHECRs 

In this subsection, we present the results of the numeri- 
cal simulations for the angular images of UHECRs in a va- 
riety of conditions concerning strength and correlation length 
of the EGMF and the limiting magnitudes of galaxies. We 
show the angular images for two energy range (10 19,6 - 10 20 ) 
eV and above 10 20 eV, about at which the GZK cutoff of the 
energy spectrum is predicted. For comparison, we also show 
the AGASA events with the angular images. It is noted that 
the current AGASA data set contains 49 events with energies 
of (10 19 ' 6 - 10 20 ) eV and 8 events above 10 20 eV. 

The calculated angular images of UHECRs originating the 
galaxies more luminous than M\ lm = -14.0 (upper), -18.0 (mid- 
dle), -20.0 (bottom) for number (left) and luminosity (right) 
weighted sources in the case of (B,l c ) = (1,1) are shown in 
Figure 5 for energy range (10 19 6 - 10 20 ) eV. The AGASA 49 
events in this energy range are also shown as circles of radius 
proportional to their energies. Spatial structure of the LSC is 
reflected by prominent high intensity region running north and 
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south at a ~ 180° forMi im = -14.0 and -18.0. Since the lumi- 
nous galaxies which distribute outward than the faint ones make 
a substantial contribution in the case of luminosity weighted 
sources, the region of high intensity becomes to be less obvious 
in the case of Mn m = -14.0. In other cases, there is no signif- 
icant difference simply due to small variance of absolute mag- 
nitudes of the individual galaxy. We can see that the angular 
image becomes to be isotropic and show good correlation with 
AGASA events as restricting sources to more luminous galax- 
ies. Such correlation of the AGASA events with angular image 
of UHECRs is also pointed out in the case of galaxies with huge 
infrared luminosity (Smialkowski et al. 2002). 

Even for Mu m = -20.0, there are high intensity regions around 
(a, S) ~ (180°, 10°) which do not correlate with the AGASA 
events. However, if we select sources, which contribute to 
the currently observed cosmic ray flux, from our ORS sample, 
these high intensity regions can be eliminated. An example in 
the case of Mr, m = -20.5 is shown in Figure 21 as we explain 
later. 

Figure 6 shows the predicted angular image of UHECRs 
for (B,l c ) =(1,1) (upper), (1,10) (middle), (10,1) (bottom) and 
M lim = -14.0 (left), -20.0 (right) with energies of (10 196 - 10 20 ) 
eV. The angular image is distorted for (B,l c ) =(1,10), (10,1), and 
dependence on the correlation length is relatively weak because 
the deflection angle of UHECRs is proportional to B ■ l l J 2 . As 
we will show later, this distortion of the angular image is crucial 
for the small scale clustering. 

In Figure 7, we show the angular image for (B, l c ) = (1, 1) (up- 
per), (10, 1) (middle), (100, 1) (bottom) and M lim =-14.0 (left), 
-20.0 (right) with the energy range (10 20 - 10 20 3 ) eV. (High- 
est energy of the events observed by the AGASA is <~ 10 20 3 
eV.) Because protons with energies above 10 20 eV must origi- 
nate within about 50 Mpc from us due to photopion production, 
only nearby galaxies are able to make a substantial contribution 
to the angular image so that the spatial structure of the LSC is 
clearly visible. There is no correlation for M\\ m = -14.0 between 
the AGASA events and the angular image of UHECRs. For 
M\ lm = -20.0, arrival distribution is relatively isotropic, but still 
has no correlation with the AGASA events. As we conclude 
later, the AGASA events above 10 20 eV might be of different 
origin, if they are confirmed. 
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FIG. 8. — The energy spectra in the case of B = 1 (upper), 10 (middle), 100 
(bottom) nG and M lim = -20 for number weighted sources. They are normal- 
ized arbitrarily, but with the same factors for all the values of the correlation 
length. 



When the correlation length of the EGMF is small, protons 
can not be deflected effectively and thus its path length and 
equivalently energy loss is smaller than that for longer corre- 
lation length. This can also be seen in Figure 4 and Figure 8. 
Especially, effects of changing the correlation length are large 
at£ < 10 20 eV forB = 10 nG and at E ~ 10 20 eV forB = 100 nG, 
where gyroradii of UHECRs are about 1 Mpc and thus deflec- 
tion is strong enough for the diffusion approximation to become 
applicable in the case of l c = 10 and 40 Mpc. 

In order to determine normalization of the energy spectrum 
and quantify statistical significance of mean deviation between 
the fitted energy spectrum and the observed one, we introduce 
Xes in the similar way as the two point correlation function ex- 
plained in the previous section. This is defined by 



Xes = ■ 



1 



(dN/dEjEQ-dN^/dEjEdY 
^ a 2 

£i=10 195 1 
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3.2. Energy Spectra of UHECRs 

We present in this subsection the results of the numerical 
simulations for the energy spectra of UHECRs. Figure 8 shows 
the energy spectrum in the case of B = 1 (upper), 10 (mid- 
dle), 100 (bottom) nG and M\ lm = -20.0 for number weighted 
sources. They are normalized arbitrary, but with the same fac- 
tors for all the values of the correlation length. The substantial 
deflection in the EGMF results in a significant increase of the 
path length, then leading to a time delay exceeding the age of 
the universe or a decrease of energy below E = 10 196 eV due 
to pair production for protons injected at sufficiently large dis- 
tance, as easily seen in Figure 4. Thus, the stronger EGMF 
make the GZK cutoff around 10 20 eV less sharp or eliminate it 
in Figure 8. 



where dN / dE \, s (Ei) is the energy spectrum observed by the 
AGASA at E = E[, and o\ is lcr error at this energy bin. 

The origin of cosmic rays with energies below 10 19,5 eV is 
also one of the major open questions in astro-particle physics. 
Berezinsky et al.(2002) showed that predicted UHECRs flux 
from GRBs fall short of the observed flux below 10 19 5 eV as- 
suming an injection spectrum E~ 2 , which is advocated in Vi- 
etri (1995) and Waxman (1995), and also used in this study. 
However, there are another possible extra-galactic and galac- 
tic UHECR production sites, such as AGNs and young neu- 
tron stars (Blasi et al. 2000). Maximum energy of cosmic 
ray achieved by them may be lower than 10 19 5 eV. In this 
case, these galactic and/or extra-galactic components may sub- 
stantially contribute to the cosmic ray flux below 10 19,5 eV. 
Throughout the paper, we assume that these components bridge 
the gap between the observed flux and the predicted one of 
Berezinsky et al.(2002), and restrict ourselves to the energy 
spectrum only above 10 19 5 eV. 
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FIG. 5. — Angular probability distribution of UHECRs with energies of (10 19 6 - 10 20 °) eV originating the galaxies more luminous than My tm = -14.0 (upper), 
-18.0 (middle), -20.0 (bottom) for number (left) and luminosity (right) weighted sources in the case of (B, l c ) = (1,1). Resolution of the image is set to be 1° X 1°. 
The AGASA 49 events in the corresponding energy range are also shown as circles of radius proportional to their energies. 
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FIG. 6. — Same as Figure 5, but for (B,k) = (1,1) (upper), (1, 10) (middle), (10, 1) (bottom) andMi im =-14.0 (left), -20.0 (right) in the case of number weighted 
sources. 
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FIG. 9. — The energy spectra in the case of B = 1 (upper), 10 (middle), 100 
(bottom) nG and M\- Im = — 14 (left), -20.0 (right) for number weighted sources. 
They are normalized so as to minimize xes. which is shown in this figure 
for l c - 1 Mpc. The energy spectrum observed by the AGASA is also shown 
(Hayashida et al. 2000). 

Figure 9 shows the energy spectrum normalized so as to 
minimize xes in the case of B = 1 (upper), 10 (middle), 100 
(bottom) nG and M Um = -14.0 (left), -20.0 (right) for number 
weighted sources. We also show the significance values xes m 
this figure. 
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FIG. 10. — Same as Figure 9, but for luminosity weighted sources. 

The predicted energy spectrum for (B,Mn m ) = (1,-14.0) pro- 
vide good fit to the AGASA data, including the extension 
above 10 20 eV. However, arrival distribution is significantly 
anisotropic and inconsistent with the AGASA (see the next sub- 
section). In the case of (B 7 M\i m ) = (1,-20.0), the energy spec- 
trum can be well fitted with the AGASA data below E = 10 20 
eV, but the GZK cutoff is predicted above E = 10 20 eV due to the 
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energy loss by photopion production in contrast to the AGASA 
data. This is because luminous galaxies near the earth are less 
than dwarf galaxies (see Figure 2). Although the energy spec- 
trum extends beyond 10 20 eV for the stronger EGMF, it does not 
provide good fit with the AGASA data at E < 10 20 eV, where 
the current data have adequate statistics. Consequently, \es f° r 
the stronger EGMF is larger than that for B = 1 nG as easily 
seen in Figure 9. 

In the same way as Figure 9, we show the energy spectra 
for luminosity weighted sources in Figure 10. Since contri- 
butions from the dwarf galaxies are suppressed in this case, 
the cutoff of the energy spectra for M\ lm = -14.0 is predicted. 
For M]j m = -20.0, the results are almost the same as number 
weighted ones. 

3.3. Statistics on the Arrival Directions of UHECRs 

In this subsection, we show the results of statistics on the 
arrival directions of UHECRs. Figure 1 1 shows the predicted 
first harmonics for number weighted sources as a function of 
limiting magnitude. The error bars represent the statistical er- 
ror due to the finite number of simulated events, whose number 
is set to be that observed by the AGASA in the energy range of 
(10 19 6 -10 20 3 ) eV. Arrival directions of UHECRs are restricted 
in the range -10° < 8 < 80° in order to compare our results 
with the AGASA data. The solid lines indicate the harmonic 
amplitude obtained by the AGASA data. The shaded region 
is expected from the statistical fluctuation of isotropic source 
distribution with the chance probability larger than 10%. 
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FIG. 11. — Amplitude of the first harmonics as a function of the limiting 
magnitude of galaxies in the case of number weighted sources. The number of 
simulated events is set to be 57 in the energy range of (10"- s - 10 20,3 ). We dial 
sets of events 1000 times from the calculated probability density distribution 
in order to obtain averages and variances of the amplitude. The solid lines rep- 
resent the AGASA data in the corresponding energy range. The shaded region 
is expected from the statistical fluctuation of isotropic source distribution with 
the chance probability larger than 10%. 

Clearly visible in Figure 1 1 is that the anisotropy is least for 
Mim ~ -20.5, which is caused by broadened distribution of lu- 
minous galaxies in the LSC. The amplitude at Mi; m = -20.5 is 
consistent with isotropic source distribution within 1 a level 
for all the (B 7 l c ). Compared with the AGASA data, a signif- 
icant anisotropy can be seen for Mr, m > -18.0, although the 



energy spectra in these cases extend above 10 20 eV as shown 
in the previous subsection. Increase of the amplitude around 
M]im <~ -21 .0 is mainly due to two giant galaxies existing in the 
almost same direction (a ~ 185°, 6 ~ 14.7°) in the vicinity of 
our Galaxy (at ~ 30 Mpc). 

Comparing the results of different strength of the EGMF, we 
find that the anisotropy in the cases of 10, 100 nG is larger than 
that of 1 nG below Mr, m ~ -18.0 contrary to expectation. We 
note again that stronger EGMF results in shorter cosmic ray 
horizon because of time delay and energy loss due to the sub- 
stantial deflection in the EGMF (See Figure 4). That is, spatial 
structure of nearby galaxies is strongly reflected in the large 
scale anisotropy for the strong EGMF. Above Mu m ~ -18.0, 
since luminous galaxies distribute outward, the anisotropy in 
the cases of the strong EGMF is smaller as expected. Also, 
increase of correlation length results in more effective deflec- 
tions, and then leads to smaller anisotropy in the case of 1, 10 
nG. However, this dependence is opposite for 100 nG, which 
is caused by a reason somewhat similar to that of slightly large 
anisotropy for the strong EGMF below Mr, m ~ -18.0 explained 
above. Gyroradius of proton with energy 10 200 eV in the 
EGMF of 100 nG is 1 Mpc, which is smaller than l c = 10,40 
Mpc. In this case, UHECR deflection is strong enough for the 
diffusion approximation to become applicable. Thus, the cos- 
mic ray horizon is extremely shorter than that for l c = 1 Mpc, so 
that the anisotropy arise from nearby galaxies in the LSC. 

Figure 12 shows the predicted second harmonics for number 
weighted sources in the same manner as the first harmonics. 
Dependence on the limiting magnitude is weaker than that of 
the first harmonics. There is no parameter set (B, l c ,Mn m ) which 
is consistent with isotropic source distribution. We also calcu- 
lated the first and second harmonics for luminosity weighted 
sources, and found that the anisotropy is slightly less than that 
of number weighted sources forMi; m > —18.0 because of small 
contributions from the dwarf galaxies, which concentrate in the 
LSC. ForMi im < -18.0, the results are almost the same as num- 
ber weighted sources. 
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FIG. 12. — Same as Figure 11, but amplitude of the second harmonics. 

The two point correlation functions are shown for number 
weighted sources in the case of Mr, m = -14.0 (left panels) and 
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M\\ m = -20.5 (right panels) in Figure 13. The histograms repre- 
sent the AGASA data, which show the statistically significant 
correlation at the smallest angle scale (~ 1°). xio defined in the 
previous section is also shown. 
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FIG. 13. — The two point correlation function for number weighted 
sources in the case of M lim = -14.0 (left panels), -20.5 (right panels) and 
(B,l c ) = (1, 1),(1, 10), (10, 1),(100, 1) in descending order. The number of sim- 
ulated events is set to be 57 with energies of (10 19 6 - 10 20 ' 3 ) eV. The his- 
tograms represent the AGASA data in this energy range, xio is also shown. 

At larger angles (~ 5°), there are relatively strong corre- 
lation for Mum = -14.0, which is not consistent with the ob- 
served isotropic distribution, whereas arrival distributions for 
Mum = -20.5 are sufficiently isotropic. Among the cases of 
Mum = -20.5, noticeable feature is coexistence of a strong cor- 
relation at small angles and sufficiently isotropic distribution 
at larger scale for (B,l c ) = (1, 1). This is a consequence of the 
sharp peak of the UHECR angular image in the case of (B, l c ) = 
(1,1) (see Figure 5). This strong correlation at small angles is 
reduced or eliminated for longer correlation length or stronger 
EGMF Due to small number of observed events, the statistics 
of the two point correlation function is limited. For this rea- 
son, xio for (B,l c ) = (1,1) is equal to that for (B,l c ) = (1, 10). 
However, if the future experiments like the Pierre Auger ar- 
ray (Capelle et al. 1998) confirm the small scale clustering of 
UHECR arrival directions with more event number, the case of 
(B,l c ) = (1, 1) will be favored with better significance. 

Let us summarize the results of the three quantities (energy 
spectrum, harmonic amplitude, two point correlation function). 
First, the energy spectrum for B = 1 nG provide better fit to 
the AGASA data than that for strong EGMF in terms of \es- 
Among this case, the extension of the energy spectrum is also 
reproduced forMvj m ~ -14.0, although there is significant large 
scale anisotropy contrary to the AGASA data in this case. In 
the case of My lm <~ -20.0, the GZK cutoff is predicted. 

Second, the amplitude of the first harmonics is consistent 
with the AGASA data at M\ im ~ -20.5 irrespective of (B,l c ). 
The amplitude of the second harmonics is not consistent with 
the AGASA for any (B, l c ,Mi im ). 

Finally, the two point correlation function for (B,l c ,Mn m ) = 
(1 , 1 ,-20.5) indicate coexistence of a strong correlation at small 



angles (~ 1°) and sufficiently isotropic distribution at larger 
scale. In sum, the source model of (B,l Cl Mu m ) = (1, 1,-20.5) 
seems to reproduce the observation better than another parame- 
ter sets. The amplitude of the second harmonics and the exten- 
sion of the energy spectrum can not be explained in this source 
model. 



3.4. Dependence on the Number of UHECR Sources 

Up to present, we have not explained the large scale isotropy 
quantified by the amplitude of the second harmonics and the 
small scale clustering except for (B,l c ) = (1, 1), as well as the 
extension of the cosmic ray spectrum above 10 20 eV. Since 
we assumed that all the galaxies more luminous than a given 
limiting magnitude contribute to the cosmic ray flux, the num- 
ber of sources is maximal and the fluctuations around the as- 
sumed non-isotropic distribution is minimal, which make the 
anisotropy most visible. Furthermore, each source contribute at 
most one event in this case, it is difficult to obtain clusters which 
are likely to reflect the point-like sources. Given this situation, 
we select some galaxies from our ORS sample, and investigate 
dependence of the results on the number of selected galaxies. 
It is also expected that we can know the effect of fluctuation 
of source number within the GZK sphere on the energy spec- 
trum above 10 20 eV. As an example, this selection corresponds 
to restricting the sources of UHECRs to host galaxies of GRBs, 
which contribute to the currently observed cosmic ray flux. 

We assumed that contributions to arrival distribution from 
sources outside 8000 km s" 1 from us, within which the ORS 
is sampling the universe, is completely isotropic one. However, 
since we select only some galaxies in the following, point-like 
nature of sources has to be taken into account. Thus, we dis- 
tribute galaxies homogeneously in this region, and including 
them, we select some galaxies from our ORS sample. We re- 
strict ourselves to the limiting magnitude of Mn m = -20.5 and 
number weighted sources, which produce the most isotropic ar- 
rival distribution of UHECRs. It is noted that the number den- 
sity of galaxies more luminous than Mr, m = -20.5 in our ORS 
sample is ~ 10 -43 Mpc" 3 . 

In Figure 14, we show the amplitude of the first harmonics 
as a function of the number fraction (NF) of selected sources 
to all the ORS galaxies more luminous than M\ lm = -20.5 for 
all the (B,l c ). For each number fraction, we plot the average 
over all trial of source selection and realization from the simu- 
lated probability distribution with two error bars. The smaller is 
the statistical error due to the finite number of observed events, 
while the larger is both the statistical error and the cosmic vari- 
ance, that is, variation between different selections of sources 
from our ORS sample. In order to obtain the average and vari- 
ance, we dial the simulated sets of events 100 times from prob- 
ability distribution predicted by a specific source distribution, 
and the sets of sources 30 and 10 times from our ORS sample 
for the number fraction of < 0.1 and > 0.1 respectively in the 
case of B = 1, 10 nG, but 30 times for all the number fractions 
in B = 100 nG. We note that the number of sources, which con- 
tributes to the cosmic ray flux, for each set of (B, l c ) differ from 
the others, even for the same number fraction, because of the 
different range of UHECRs (See Figure 4). For example, there 
is no source within 100 Mpc for the number fraction of < 10~ 2 5 
(for the number density of < 10" 7 Mpc" 3 ) in the case of B = 100 
nG. 
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FIG. 14. — Amplitude of the first harmonics as a function of the num- 
ber fraction of UHECR sources relative to the ORS galaxies luminous than 
M]j m = -20.5. UHECR sources are selected from the ORS sample of M\\ m < 
-20.5. The number of simulated events is set to be 57 in the energy range 
of (10 19 6 - 10 20 ' 3 ). The points indicate the realization average, and the left 
and right error bars represent the statistical and total error, respectively. The 
solid line represent the AGASA data in the corresponding energy range. The 
shaded region is expected from the statistical fluctuation of isotropic source 
distribution with the chance probability larger than 10%. 

As a general trend, anisotropy increases with decreasing the 
number fraction, because the arrival directions of UHECRs 
clumps to that of sources which is small in number. For large 
number fraction (~ 1), the amplitudes approach to the values 
in the case that all galaxies more luminous than M\{ m = -20.5 
contribute to the observed cosmic ray flux, as they should. 
As clearly visible in Figure 4, UHECRs can propagate the 
longest distance for (B,l c ) = (1, 1), in which the deflection an- 
gle, and then the effective path length, is smallest among for 
another parameters. Accordingly, distant sources make a sub- 
stantial contribution to the arrival distribution of UHECRs, un- 
like in the case of stronger EGMF This is reflected in the de- 
crease of the amplitude around the number fraction of 10 -1 7 
for (B, l c ) = (1, 1). Except for B = 100 nG, arrival distribution is 
sufficiently isotropic compared to isotropic source distribution 
for almost all the number fractions. 

The amplitude of the second harmonics is also shown in Fig- 
ure 15. Dependence on the number fraction is roughly same as 
that of the first harmonics. For B = 1 and 10 nG, the predicted 
amplitude is sufficiently isotropic for the number fraction of 
(10~ 30 - 10" 1 5 ). On the other hand, it appears to be inconsis- 
tent with the observation for any number fraction in the case of 
B = 100 nG. 

Figure 16 shows xio of the two point correlation function 
as a function of the number fraction. Increase of deviations 
between models and data with decreasing the number fraction 
is attributed to the clump of UHECRs at directions of their 
sources, similarly to the case of the harmonic analysis. Es- 
pecially, this dependence is quite noticeable for (B,l c ) = (1, 1), 
where UHECRs hardly deflect so that correlation at small an- 
gle is too strong. At larger number fraction (> 10~ 2 °), it can be 
seen that the model of (B, l c ) = (1, 1) provide relatively better fit 
to the observation than another models. 
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FIG. 15. — Same as Figure 14, but amplitude of the second harmonics. 

In Figure 17, we show the two point correlation function pre- 
dicted by a specific source scenario for (B,l c ) = (1,1), (1, 10) 
and (10,1) with the number fraction of <~ 10~ 1,7 . xio is alS0 
shown. Each source distribution is selected so that it predicts 
the smallest value of xio for respective set of (B,l c ). Only for 
(B,l c ) = (1,1), another condition is also imposed that the har- 
monic amplitude predicted by the source model is consistent 
with the AGASA observation within la level. The correlation 
of events at small angle scale and sufficient isotropic distribu- 
tion at larger scale are well reproduced for (B, l c ) = (1, 1). There 
is no significant correlation for (B,l c ) = (1, 10) and (10, 1), be- 
cause of the larger deflection angle of UHECRs. In the cases 
of strong EGMF or longer correlation length, it appears to be 
difficult to reproduce the strong anisotropy only at the smallest 
angle scale, even for the case of minimal value of Xio- 
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FIG. 16. — xio as a function of the number fraction of UHECR sources 
relative to the ORS galaxies luminous than M\ lm = -20.5. UHECR sources 
are selected from the ORS sample of M\\ m < -20.5. The errorbars represent 
cosmic variance due to different realizations of the source selection. 
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FIG. 17. — The two point correlation function predicted by a specific source 
scenario in the case that the number fraction (NF) ~ 10 _L7 of the ORS galaxies 
more luminous than M\ lm = -20.5 is selected as UHECR sources. The number 
of simulated events is set to be 57 with energies of (10 19 - 6 - 10 203 ) eV. The 
histograms represent the AGASA data in this energy range, xio is also shown. 

Considering the calculated results of the first and second har- 
monics and the small scale anisotropy all together, the model of 
(B,/ c ) = (1,1) and the number fraction of ~ 10~ L7 (the source 
number density of ~ 10~ 6 Mpc" 3 ) seems to reproduce the ob- 
servations better than another parameter sets. However, due to 
the small number of observed events, the statistics for the two 
point correlation function is limited. For this reason, xio does 
not differ very much from each other. The future experiments 
like the Pierre Auger array (Capelle et al. 1998) will decrease 
statistical uncertainty and provide more strong constraints to 
the model predictions. 



Next we discuss dependence of the energy spectrum on the 
number fraction. In Figure 18, we show the predicted energy 
spectra in the case of (B,/ c ) = (1,1) for the number fraction 
lO" 07 , 1(T 17 , 1(T 2 7 in descending order. They are normal- 
ized so as to minimize \es- The shaded regions represent 
1 a error due to the source selection from our ORS galaxies 
more luminous than M\ im = -20.5. As is evident from this fig- 
ure, decreasing the number fraction increases fluctuation of the 
predicted energy spectrum, especially above 10 20 eV. For the 
number fraction of 10~ 2 7 (the source number density of ~ 10~ 7 
Mpc" 3 ), since there are no source within the GZK sphere, the 
GZK cutoff is clearly visible. The predicted energy spectrum 
is roughly consistent with the AGASA observation within 1 a 
level in the case of the number fraction of 10~ 1,7 , where the 
observed isotropy and clusters are also reproduced well. 

However, this result should be interpreted with care. Since 
the source number density is ~ 10~ 6 Mpc -3 in the case of the 
number fraction of 10~ 1,7 , mean number of sources is ~ 0.5 
within the GZK sphere. Therefore, the AGASA 8 events above 
10 20 eV, which do not constitute the clustered events with each 
other, must originate from at most a few sources. We found a 
source model which explain the extension of the energy spec- 
trum. In this source model, there are two sources within the 
region covered by the ORS (~ 107 Mpc), one at ~ 20 Mpc, the 
other at ~ 105 Mpc. In Figure 19, we show the energy spec- 
trum predicted by this source model. The dotted line repre- 
sent contribution from the nearest source. The predicted energy 
spectrum provide good fit to the observed one including the ex- 
tension beyond 10 20 eV. However, almost all the events above 
10 20 are generated by the nearest source as evident from Fig- 
ure 19. These events are strongly concentrated in a single cen- 
ter, which is inconsistent with the AGASA 8 events. It seems to 
be difficult to reproduce the AGASA observation above 10 19 6 
eV (including above 10 20 eV) by a single scenario of UHECR 
origin, even for (B,l c ,M\\ m ,NF) = (1, 1,-20.5, 10" L7 ). 
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FIG. 18. — The energy spectra predicted by sources selected from the ORS 
galaxies more luminous thanMi; m =-20.5 in the case of (B, l c ) = (1, 1). NF rep- 
resent the number fraction of selected UHECR sources to all the ORS galaxies 
Wim < -20.5). They are normalized so as to minimize Xes- The shaded re- 
gions represent 1 a error due to the source selection from our ORS sample. 
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FIG. 19. — The energy spectrum predicted by a specific source scenario 
in the case that the number fraction (NF) ~ 10 - of the ORS galaxies more 
luminous than M[- Im = —20.5 is selected as UHECR sources. The dotted line 
represent contribution from the nearest source (at ~ 20 Mpc) in the source 
model. The energy spectrum is normalized so as to minimize xes- The energy 
spectrum observed by the AGASA is also shown. 
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FIG. 20. — Same as Figure 17, but for the number fraction of ~ 10 '. 

Recently, the HiRes Collaboration (Wilkinson et al. 1999) 
reported the cosmic ray flux from 2 x 10 17 eV to over 10 20 eV 
including the GZK cutoff (Abu-Zayyad et al. 2002). In our nu- 
merical simulations, the sharp GZK cutoff of the energy spec- 
trum is predicted for the number fraction of ~ 10" 2 7 as shown 
above. The small scale correlation of events in this case is too 
strong for (B, l c ) = (1, 1) as is evident from Figure 20, where the 
two point correlation function is shown in the same way as Fig- 
ure 17, but with the number fraction of <~ 10" 2 1 . This strong 
correlation is reduced for (B,l c ) = (1, 10) and (10, 1), but still 
inconsistent with the AGASA data. 




FIG. 21. — Arrival directions of UHECRs with energies of (10 1 ' 6 - 
10 2 °3) £ y m me e q Ua t or j a i coordinate predicted by the source model of 
(B,l c ,M hm ,NF) = (1, 1,-20.5, Kr 17 ) of Figure 17. For a comparison, the 
AGASA 57 events in this energy range are also shown. The number of sim- 
ulated events is also set to be 57. The arrival directions are restricted in the 
range —10° < S < 80° in order to compare our results with the AGASA data. 
r m represent the m-th harmonic amplitude. 

Finally we show in Figure 21 realizations of 57 arrival direc- 
tions of UHECRs with energies of (10 19 6 - 10 20 3 ) eV, predicted 
by the source model of (B,l c ,M Um ,NF) = (1, 1,-20.5, 10" 1 7 ) 
of Figure 17. The AGASA 57 events in this energy range 



are also shown. The arrival directions are restricted in the 
range -10° < <5 < 80° in order to compare our results with 
the AGASA data. The arrival distribution of UHECRs seems 
to be isotropic on a large scale with clusters of events mean- 
ing the small scale anisotropy. The amplitudes of the first and 
second harmonics are consistent with the AGASA, but the ex- 
tension of the spectrum is not. (Mean number of simulated 
events above 10 20 eV is only 1.) In this example, C6 doublet 
at (a, <5) ~ (210°, 35°) is well reproduced. Also, high intensity 
regions around (a, S) ~ (180°, 10°) mentioned in Section 3.1. 
are eliminated under the favor of the source selection. 

4. DISCUSSION 
4. 1 . Constraint on UHECR source model 

We presented the results of numerical simulations on the 
propagation of UHE protons with energies of (10 19 5 - 10 22 ) eV 
injected by a discrete distribution of galaxies using the ORS 
data sample, which is properly corrected for the selection effect 
and absence of galaxies in the zone of avoidance (\b\ < 20°). 
We adopted strength of the EGMF to be not only 1 nG but also 
10,100 nG, in order to explore the possibility of reproducing the 
extension of the energy spectrum above 10 20 eV. We calculated 
the three observable quantities, cosmic ray spectrum, harmonic 
amplitude, and two point correlation function. With these quan- 
tities, we made comparisons between the model predictions and 
the existing AGASA data. 

At first, we explored the source model consistent with the 
current observations, as a function of the limiting magnitudes 
of galaxies. The energy spectrum for the strong EGMF shows 
the extension beyond 10 20 eV. However, it can not be well fit- 
ted with the observed energy spectrum, especially at E < 10 20 
eV where the current data have adequate statistics. On the other 
hand, the energy spectrum for B = 1 nG provide good fit with the 
observation below E = 10 20 eV, and statistical significance of 
deviation xes is smaller than that for the strong EGMF. Among 
this case, the extension of the energy spectrum is also repro- 
duced forMii m ~ -14.0, although there is significant large scale 
anisotropy contrary to the AGASA data in this case. In the case 
of Mi im ~ -20.0, the GZK cutoff is predicted. 

We found that the angular image became to be isotropic 
and showed good correlation with AGASA events as restrict- 
ing sources to more luminous galaxies. As a consequence, 
galaxies more luminous than -20.5 mag produce the angular 
image of UHECR which is most isotropic roughly irrespective 
of strength and correlation length of the EGMF. However, it is 
not isotropic enough to be consistent with the observed ampli- 
tude of the second harmonics, even forMi; m = -20.5. We further 
found that the small scale anisotropy is well reproduced in the 
case of (B,l c ) = (1,1) because of the small deflection angle of 
UHECRs. 

Next, in order to obtain sufficiently isotropic arrival distribu- 
tion of UHECRs and small scale clustering even for the case of 
other than (B,l c ) = (1,1), we randomly selected sources from 
our ORS sample, and investigate the dependence of the re- 
sults on their number. We found that the arrival distribution 
became to be sufficiently isotropic in the case that the number 
fraction (10~ 3 - 10~ L5 ) of the ORS galaxies more luminous 
than M\{ m = -20.5 is selected as UHECR sources for B = 1, 10 
nG. Among this case, the small scale clustering can also be re- 
produced for (B,l c ,Miim,NF) = (1, 1,-20.5, 10" 1 7 ) better than 
the cases of another parameter sets. We furthermore found that 
the extension of the energy spectrum is explained in this case 
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within roughly 1 a error due to source selection from our ORS 
sample because of the large fluctuation above 10 20 eV. 

In sum, we showed that the three observable quantities in- 
cluding the GZK cutoff of the energy spectrum can be repro- 
duced in the case that the number fraction ~ 10~ 1,7 of the ORS 
galaxies more luminous than -20.5 mag is selected as UHECR 
sources. In terms of the source number density, this constraint 
corresponds to ~ 10" 6 Mpc" 3 . 

4.2. Implication for sources of UHECR above 10 20 eV 

As mentioned in the previous subsection, we showed that 
the three observable quantities including the GZK cutoff of the 
energy spectrum can be reproduced in the case that the num- 
ber fraction ~ 10 -1 7 of the ORS galaxies more luminous than 
-20.5 mag is selected as UHECR sources. However, we should 
interpret this result with care. Since the AGASA 8 events above 
10 20 eV do not constitute the clustered events with each other, 
the number of their sources within the GZK sphere must be 
much larger than the event number (8). On the other hand, the 
number of sources within the GZK sphere is at most only 1 
in the case of the number fraction of 10~ L7 (the source num- 
ber density of 10~ 6 Mpc" 3 ). It is suspected that cosmic rays 
above 10 20 eV (hereafter EHECRs; extremely-high energy cos- 
mic rays) are strongly concentrated in a single center, even if 
the extension of the energy spectrum is explained. Therefore, 
large fraction of EHECRs observed by the AGASA might orig- 
inate from sources other than ones of UHECR with energies of 
(10 19 6 -10 20 ) eV. The number density of EHECR source must 
be much larger than 10~ 6 Mpc" 3 . 

We can put forward two candidates as such EHECR produc- 
tion site. It is suggested that iron ions from the surfaces of 
young strongly magnetized neutron stars may be accelerated 
to 10 20 eV through relativistic MHD winds (Blasi et al. 2000). 
However, it seems to be difficult to obtain sufficiently isotropic 
arrival directions of EHECRs from sources localized in the 
Galactic plane by deflection due to the Galactic magnetic field 
(O'Neill etal. 2001). 

The decay of some supermassive particles, which could be 
produced from TDs, or be certain MSRPs, are also considered 
as probable EHECRs origin. Decay rate of the supermassive 
particles required to explain the observed EHECR flux are es- 
timated as — 10 35 Mpc" 3 yr" 1 (Bhattacharjee and Sigl 2000). 
Multiplying by the current EHECR observing time (~ 10 yr), 
we obtain the number density of <~ 10 36 Mpc" 3 , which is suf- 
ficiently large to reproduce not clustered AGASA 8 EHECRs. 
As for the large scale anisotropy, distribution of TDs in the uni- 
verse would be homogeneous. MSRPs are expected to be clus- 
tered in galactic halos, and thus EHECR flux will be dominated 
by contribution from our own Galactic Halo. Accordingly, ar- 
rival distribution would be isotropic enough to be consistent 
with the current AGASA observation in both cases. Further- 
more, these scenarios are not constrained by contribution to the 
cosmic ray flux at lower energy which must be smaller than 
observed one, since hard injection spectrum is generally pre- 
dicted. For the reasons stated above, the top-down scenarios 
may be favored to explain large fraction of the 8 EHECRs ob- 
served by the AGASA. 

There is a problem if all the EHECRs observed by the 
AGASA originate from sources other than those of UHECRs. 
If we consider that sources of EHECRs and UHECRs precisely 
differ from each other, there would be no clustered event set 
which includes both EHECRs and UHECRs, because clustered 



events are expected to reflect a single point source. However, 2 
EHECRs out of 8 observed by the AGASA constitute the clus- 
tered events with UHECRs with energies of (10 19 6 - 10 20 ) eV 
(CI and C3 in Hayashida et al.(2000)). We have to interpret 
that the two clusters happen to be produced by chance. Chance 
probability of more than two accidental clusterings within 2.5° 
between 49 UHECRs and 8 EHECRs is quite low (~ 4%). The 
AGASA observation seems to deviate from our prediction with 
~ 2a confidence level. 

However, about 1 EHECR out of 8 may originate in the 
bottom-up scenarios. Indeed, 57 cosmic rays above 10 19,6 eV 
predicted by the source model of Figure 21 include an EHECR. 
In this case, since the source of this single event will be located 
within the GZK sphere, it is likely that this event constitutes a 
clustered event set with UHECRs. Chance probability of more 
than one accidental clustering rise up to ~ 28%, which is high 
enough for us to interpret one of the two doublets to be pro- 
duced merely by chance. Thus, our scenario is consistent with 
the current AGASA observation within 1 a level. 

The event correlation between EHECRs and UHECRs would 
provide an important test of our conclusion with sufficient 
amount of data from future experiments. For example, when the 
event number of UHECRs above 10 19 ' 6 increases to 10 times as 
many as that of the current AGASA observation (57), the num- 
ber of EHECRs which lie within 2.5° from UHECRs by chance 
becomes 18.6 ±4.3 (la). If future experiments observe event 
number of such EHECRs within this range, our scenario will 
be supported. Here we note that the source number within the 
GZK sphere is at most 1 in the case of NF=10 _1 1 . Therefore the 
number of event clustering between an EHECR and an UHECR 
which can be explained by bottom-up scenarios is at most 1 . 

Another test of our conclusion is distinguishing the chem- 
ical composition of primary cosmic rays which initiate air 
showers in the earth atmosphere. In the bottom-up scenarios, 
protons or nuclei are primary particles. On the other hand, 
photons are considered as primary in the top-down scenar- 
ios (Bhattacharjee and Sigl 2000). Information on the chem- 
ical composition is mainly provided by the depth of shower 
maximum for fluorescence observation of the air shower. 
There are a number of projects utilizing the fluorescence tech- 
nique, such as HiRes (Wilkinson et al. 1999), South and North 
Auger (Capelle et al. 1998), OWL (Cline and Stecker 2000) 
and EUSO (Benson and Linsley 1992). These experiment 
would provide information on the chemical composition and 
a test of our conclusion. 

4.3. Comparison with the HiRes cosmic ray spectrum 

The HiRes Collaboration (Wilkinson et al. 1999) reported 
the cosmic ray flux from 2 x 10 17 eV to over 10 20 eV includ- 
ing the GZK cutoff (Abu-Zayyad et al. 2002). As discussed 
above, we showed that the three observable quantities except 
for EHECR events observed by the AGASA can be repro- 
duced in the case that the number fraction <~ 10" 1 7 of the ORS 
galaxies more luminous than -20.5 mag is selected as UHECR 
sources. In Figure 22, we show the predicted energy spectrum 
fitted to the data measured by Hires-I detector in the case of 
(B,l c ,M Umi NF) = (1, 1,-20.5, 10" 17 ). As easily seen, the fit is 
good as compared with the middle panel of Figure 18. If the 
flux measured by the HiRes Collaboration is correct and obser- 
vational features about the arrival distribution are same as the 
AGASA, our source model can explain both the arrival distribu- 
tion of UHECRs and the flux at the same time. When the large- 
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area detectors like the Pierre Auger array (Capelle et al. 1998) 
start operating, the statistical error of the observed energy spec- 
trum will decrease. This allows us to draw a conclusion about 
the origin of EHECRs or whether they exist or not. 
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FIG. 22. — The energy spectra predicted by sources selected from the 
ORS galaxies more luminous than M llm = -20.5 in the case of (B,l c ,NF) = 
(1, 1,-20.5, 10" 17 ). NF represent the number fraction of selected UHECR 
sources to all the ORS galaxies (Mij m < -20.5). They are fitted to the data of 
HiRes-I detector (squares and errorbars). The shaded regions represent 1 a 
error due to the source selection from our ORS sample. 



4.4. Implication for sources of UHECR below 10 20 eV 

In this subsection we discuss implications of our results, ob- 
tained by comparison of the model predictions with the exist- 
ing AGASA data, for the origin of UHECRs with energies of 
(10 19 6 - 10 20 ) eV. The number densities of possible UHECR 
production sites and ratios to our result (~ 10" 6 Mpc" 3 ) are tab- 
ulated in Table 1. 

Blast waves in AGN jet could accelerate particle to a few 
tens of EeV (Halzen and Zas 1997) and similarly for AGN 
cores. However, due to the interaction of UHECRs with very 
high radiation fields in and around the central engine of an 
AGN, these maxima are unlikely to be achieved. This is con- 
sistent with our result, because the number density of AGN 
are much larger than our results, and thus, arrival distribution 
may be inconsistent with the AGASA if all the AGNs con- 
tribute to the cosmic ray flux above 10 19 6 eV. Among AGNs, 
the most promising acceleration sites for UHECRs are the 
so-called hot-spots of FR-II radio galaxies (Biermann 1997; 
Norman et al. 1995). However, their number density seems 
to be too small, as compared with our result. It is unlikely 
that FR-II radio galaxies are responsible for all the AGASA 
49 events with energies of (10 19 6 - 10 20 ) eV. Another possible 
UHECRs production sites are BL Lacertae objects (BL Lacs) 
(Tinyakov and Tkachev 2001). Their number are larger than 
FR-II radio galaxies, but still insufficient to reproduce all the 
49 events. 

GRB can be shown to in principle accelerate protons up 
to 10 20 eV (Waxman 2000). The effective density of GRB 
is determined by the GRB rate (~ 2 x 10" 10 h 3 Mpc" 3 yr" 1 

(Mao and Mo 1998)) and the typical time delay r of UHECRs, with our result ~ 10" 6 Mpc 



which depends on the propagation distance from sources to the 
earth. Since we are now considering UHECRs with energies of 
(10 19,6 - 10 20 ) eV, we have to include contribution from sources 
outside the GZK sphere when calculating the typical time de- 
lay. We calculated the mean time delay from our numerical 
data assuming homogeneous GRBs distribution, and obtained 
r - 10 6 yr. This gives the GRB effective density ~ 10" 4 Mpc -3 . 
From our result, if about 1 % of GRBs produce UHECRs, ar- 
rival distribution predicted by GRBs would be consistent with 
the AGASA. 

The compact dynamo model has been proposed as 
a natural mechanism for accelerating cosmic rays in 
dead quasars (Boldt and Ghosh 1999; Levinson 2000; 
Boldt and Loewenstein 2000). Only 0.2 % of dead quasars 
are required to contribute to the comic ray flux above 10 19,6 eV. 

Large fraction of luminous infrared galaxies (LIRGs) are 
found to be interacting systems, suggesting that collisions and 
merging processes are responsible for triggering the huge IR 
light emission (Sanders et al. 1988). In such systems, favor- 
able environments for accelerating cosmic rays to 10 20 eV are 
provided by amplified magnetic fields on the scale of tens kpc 
resulting from gravitational compression, as well as high rel- 
ative velocities of galaxies and/or superwinds from multiple 
supernovae explosions (Cesarsky 1992). According to Smi- 
alkowski et al.(2002), we take the number density of such col- 
liding galaxies as that of galaxies with huge far infrared lumi- 
nosity Lfir > 1O H L0. The number density is roughly equal to 
ourresult. If all the galaxies with Lfir > lO n Z,0 are responsible 
for UHECRs, predicted arrival distribution would be consistent 
with the AGASA. Indeed, Smialkowski et al.(2002) showed 
that arrival distribution of UHECRs predicted by LIRGs is suf- 
ficiently isotropic such that statistical tests they used for large 
scale anisotropy was not conclusive for distinguishing between 
isotropic and LIRGs source distribution. Furthermore, they 
showed that the probability of occurring the clustered events 
are more than 10 times higher for LIRGs source distribution 
than that for isotopic distribution. Their results well agree with 
ours. 



4.5. Consistency with a statistical analysis of event cluster 

Our result is supported by a statistical analysis of cluster- 
ing of UHECRs developed by Dubovsky et al.(2000), which 
allow an estimate of the minimum number density of sources. 
They applied their arguments to the 14 events above 10 20 eV 
observed by several experiments. They obtained the minimum 
number density of sources ~6x 10~ 3 Mpc" 3 . Since we are 
now considering UHECRs below 10 20 eV, we have to apply the 
arguments given in their paper to the cosmic ray events below 
10 20 eV. 

Let us apply the arguments to the 49 events with energies 
of (10 196 - 10 200 ) eV observed by the AGASA. According to 
them, the number density must be larger than 



1 AT 3 
* 4fl 3 JV 2 ' 



(12) 



where 7Y tot , N c \, R are the total number of events, the number 
of events in clusters, e-folding length of cosmic ray luminos- 
ity from a single source respectively. The value of R can be 
obtained as ~ 700 Mpc from Figure 4. Taking 7Y tot = 49 and 
N c \ = 9, Eq.( 12) gives /z* ~ 10~ 6 Mpc" 3 which is consistent 
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Table 1 

Number density of possible UHECR production site 



Obiect 


Densitv 
[Mpc -3 ] 


Ratio 
to our result 
(10 -6 Mpc -3 ) 


ref 


AGN 


1 x 10 -4 


100.0 


(Loveday et al. 1992) 


FR-II radio galaxy 


3 x 1(T 8 


0.03 


(Woltjer 1990) 


BL Lac 


3 x 10 -7 


0.3 


(Woltjer 1990) 


GRB 


1 x 1(T 4 


100.0 


(Mao and Mo 1998) 


dead quasar 


5 x 10 -4 


500.0 


(Boldt and Ghosh 1999) 


colliding galaxy 


7 x 10 -7 


0.7 


(Smialkowski et al. 2002) 



4.6. Strong EGMF in the LSC 

Throughout the paper, we assumed that the statistical prop- 
erties of the EGMF are uniform. However, simple analyt- 
ical arguments based on magnetic flux freezing, and large 
scale structure simulations passively including the magnetic 
field (Kulsrud et al. 1997) demonstrate that the magnetic field 
is most likely as structured as are the baryons. The EGMF 
as strong as ~ l^iG in sheets and filaments of large scale 
galaxy distribution, such as in the LSC, are compatible with 
existing upper limits on Faraday rotation (Ryu et al. 1998; 
Blasi et al. 1999). It is suspected that the arrival distribution of 
UHECRs depends on the fields in the immediate environment 
of the observer. 

However, numerical simulations of UHECR propagation in 
inhomogeneous EGMF over cosmological distances are highly 
time-consuming. With the assumption of homogeneous EGMF, 
we have performed numerical simulation of spherically sym- 
metric propagation. Provided that the EGMF is inhomoge- 
neous, the propagation of UHECRs from a single point source 
has no longer spherical symmetry. As a result, we have to spec- 
ify the earth position in the universe. We also have to choose 
the detector (earth) size so small enough for us to accurately 
calculate the arrival directions. In this case, the number fraction 
of UHECRs arriving at the earth to injected ones is extremely 
small. This requires the number of particle to be propagated 
several orders of magnitudes higher than that used in this study, 
which takes enormous CPU time. 

Although the strong EGMF in the LSC is compatible with 
the Faraday ratation measures, the cosmic ray observations do 
not appear to support such local EGMF. Indeed, we showed in 
this study that small scale clustering can be well reproduced in 
the case that UHECRs propagate along nearly straight lines. If 
the local strong EGMF affects the arrival directions of UHE- 
CRs, small scale clustering observed by the AGASA will not 
be obtained. 

This difficulty can be seen from the studies by another au- 
thors. There are several groups who study the propagation of 
UHECR in such strong EGMF in the LSC (Sigl et al. 1999; 
Lemoine et al. 1999; Isola et al. 2002; Sigl 2002). They as- 
sume simplified source distributions which represent the LSC. 
The strong EGMF of ~ 1/iG leads to substantial deflections 
of UHECRs, which are better for explaining the observed 
isotropic arrival distribution of UHECRs. However, consis- 
tency of small scale anisotropy and also large scale isotropy 



predicted by their scenarios with the AGASA observation is 
somewhat worse than that predicted by our scenario. The strong 
EGMF in the LSC seems to be unfavored by the current cosmic 
ray observation, although it is not ruled out by the Faraday ro- 
tation measurements. 

Of course, since the event number of UHECRs is currently 
very small, we can not absolutely rule out such strong EGMF 
in the LSC. Fulture experiments would increase the significance 
level of statement about whether such strong EGMF affects the 
arrival directions of UHECRs or not. 

5. CONCLUSION 

In this section, we present our conclusions obtained from 
comprehensive study on the bottom-up scenarios. For the origin 
of EHECRs, we conclude that large fraction of the AGASA 8 
events above 10 20 eV might originate in the topdown scenarios, 
or that the cosmic ray flux measured by the HiRes experiment 
might be better. We also discussed the origin of UHECRs below 
10 20 eV through comparisons between the number density of 
astrophysical source candidates and our result (~ 10" 6 Mpc -3 ). 
At present, we can not conclude which source candidates are 
responsible for UHECRs below 10 20 eV. However, our result 
for the source number density (~ 10~ 6 Mpc -3 ) will be useful 
when searching for the ultimate UHECR source. 

Throughout the paper, we set constraints to the source mod- 
els so that their predictions reproduce the current AGASA ob- 
servation within 1 a error of the statistical fluctuation due to 
finite number of the observed events. However, if we relax the 
constraints to within 2 a error, the energy spectrum (Figure 9 
and 10) and the amplitude of the second harmonics (Figure 12) 
can be reproduced in the case of (Z?,M ]im ) ~ (1,-20.5) with- 
out the source selection. In this case, since the source number 
density is ~ 5.4 x 10" 5 Mpc -3 , there are ~ 30 sources within 
the GZK sphere. The AGASA events above 10 20 eV, which 
do not constitute the clustered events with each other, may also 
be reproduced, contrary to the case of the source number den- 
sity of ~ 10 -6 Mpc -3 . Thus, we can not absolutely rule out the 
bottom-up scenarios as the origin of EHECRs. 

The statistics for the quantities which we used is limited 
by the small number of observed events, which make dif- 
ficult to draw a decisive conclusion on the origin of UHE- 
CRs. However, the present working or development of large- 
aperture new detectors, such as HiRes (Wilkinson et al. 1999) 
and South and North Auger (Capelle et al. 1998) will con- 
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siderably decrease the statistical uncertainties. There 
are furthermore plans for space based air shower de- 
tectors such as OWL (Cline and Stecker 2000) and EUSO 
(Benson and Linsley 1992) which would detect about 1000 
events per year above 10 20 eV (Bhattacharjee and Sigl 2000). 
With the method to constrain the source model and their num- 
ber density developed in this paper, these experiments will re- 
veal the origin of UHECRs with better significance in the very 
near future. 
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